Regional Temperature Trends

Regional Surface and Bottom Temperature Trends of Northeast US

Determining rate of change and long-term net-change for surface and bottom temperatures of the Northeast US shelf regions.

library(raster)

Loading required package: sp

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::extract() masks raster::extract()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks stats::lag()
✖ dplyr::select()  masks raster::select()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(gmRi)
library(rnaturalearth)
library(scales)

Attaching package: 'scales'

The following object is masked from 'package:purrr':

    discard

The following object is masked from 'package:readr':

    col_factor

library(sf)

Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

library(gt)
library(patchwork)

Attaching package: 'patchwork'

The following object is masked from 'package:raster':

    area

library(tidyquant)

Loading required package: PerformanceAnalytics
Loading required package: xts
Loading required package: zoo

Attaching package: 'zoo'

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric


######################### Warning from 'xts' package ##########################
#                                                                             #
# The dplyr lag() function breaks how base R's lag() function is supposed to  #
# work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
# source() into this session won't work correctly.                            #
#                                                                             #
# Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
# conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
# dplyr from breaking base R's lag() function.                                #
#                                                                             #
# Code in packages is not affected. It's protected by R's namespace mechanism #
# Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
#                                                                             #
###############################################################################

Attaching package: 'xts'

The following objects are masked from 'package:dplyr':

    first, last


Attaching package: 'PerformanceAnalytics'

The following object is masked from 'package:graphics':

    legend

Loading required package: quantmod
Loading required package: TTR
Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

conflicted::conflict_prefer("select", "dplyr")

[conflicted] Will prefer dplyr::select over any other package.

conflicted::conflict_prefer("filter", "dplyr")

[conflicted] Will prefer dplyr::filter over any other package.

# Plot theme
theme_set(theme_gmri(rect = element_rect(fill = "white", color = NA), legend.position = "bottom"))

# Degree symbol
deg_c <- "\u00b0C"

# vectors for factor levels
area_levels <- c("GoM", "GB", "SNE", "MAB")
area_levels_long <- c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")

# Monthly Bottom Temperature: Source 00_Bottom_Temp_Processing.R
bpath <- cs_path("res", "Du_Pontavice_Combined_BT/RegionalTimeseries")
all_bt_monthly <- read_csv(str_c(bpath, "trawl_region_monthly_bottom_temps.csv"))

Rows: 2976 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (2): survey_area, month
dbl  (2): bot_temp, year
date (1): date

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Daily Surface Temperature: Source github.com/adamkemberling/oisst_mainstays
gb_sst  <- oisst_access_timeseries(region_family = "trawl regions", poly_name = "georges bank", box_location = "cloudstorage")
gom_sst <- oisst_access_timeseries(region_family = "trawl regions", poly_name = "gulf of maine", box_location = "cloudstorage")
sne_sst <- oisst_access_timeseries(region_family = "trawl regions", poly_name = "southern new england", box_location = "cloudstorage")
mab_sst <- oisst_access_timeseries(region_family = "trawl regions", poly_name = "mid atlantic bight", box_location = "cloudstorage")

# Combine and average monthly
all_sst_monthly <- bind_rows(
  list(
    "Gulf of Maine" = gom_sst,
    "Georges Bank" = gb_sst,
    "Southern New England" = sne_sst,
    "Mid-Atlantic Bight" = mab_sst), 
  .id = "survey_area") %>% 
  group_by(
    survey_area,
    year = lubridate::year(time),
    month = lubridate::month(time),
    month = str_pad(month, side = "left", pad = "0", width = 2)) %>% 
  summarise(sst = mean(area_wtd_sst, na.rm = T),
            .groups = "drop")

# Combine Surface and Bottom
regional_temps <- full_join(
  select(all_bt_monthly, survey_area, year, month, bot_temp),
  all_sst_monthly) %>% 
  mutate(
    date = as.Date(str_c(year, month, "15", sep = "-")),
    survey_area = fct_relevel(survey_area, area_levels_long))

Joining with `by = join_by(survey_area, year, month)`

# Add anomalies from 1991-2020
regional_clim <- regional_temps %>% 
  filter(year %in% c(1991:2020)) %>% 
  group_by(survey_area, month) %>% 
  summarise(clim_sst = mean(sst),
            clim_bt = mean(bot_temp)) %>% 
  full_join(regional_temps) %>% 
  mutate(sst_anom = sst - clim_sst,
         bt_anom = bot_temp - clim_bt)

`summarise()` has grouped output by 'survey_area'. You can override using the
`.groups` argument.
Joining with `by = join_by(survey_area, month)`

# Load polygons for aggregate regions based on trawl survey strata
regions <- read_sf(here::here("Data/raw", "nmfs_trawl_regions_collection.geojson"))


# Prepare the regions as single polygons:
gb <- filter(regions, area == "Georges Bank") %>% st_union() %>% st_as_sf()
gom <- filter(regions, area == "Gulf of Maine") %>% st_union() %>% st_as_sf()
mab <- filter(regions, area == "Mid-Atlantic Bight") %>% st_union() %>% st_as_sf()
sne <- filter(regions, area == "Southern New England") %>% st_union() %>% st_as_sf()

# month weighting so we don't have to go back to daily when averaging
seas_month_wts <- data.frame(
    "month" =  c("03", "04", "05", "09", "10", "11"),
    "n_days" = days_in_month(
      as.Date(c("2000-03-01","2000-04-01","2000-05-01",
                "2000-09-01","2000-10-01","2000-11-01"))))


# Calculate average for the season
seasonal_temps <- regional_clim %>% 
  left_join(seas_month_wts) %>% 
  mutate(
    season = ifelse(
      month %in% str_pad(c(3:5), side = "left", pad = "0", width = 2), 
      "Spring", NA),
    season = ifelse(
      month %in% str_pad(c(9:11), side = "left", pad = "0", width = 2), 
      "Fall", season),
    survey_area = factor(survey_area, levels = area_levels_long)) %>% 
  drop_na(season) %>% 
  group_by(survey_area, year, season) %>% 
  summarise(
    sst = weighted.mean(sst, w = n_days),
    clim_sst = weighted.mean(clim_sst, w = n_days),
    sst_anom = weighted.mean(sst_anom, w = n_days),
    bot_temp = weighted.mean(bot_temp, w = n_days),
    clim_bt = weighted.mean(clim_bt, w = n_days),
    bt_anom = weighted.mean(bt_anom, w = n_days),
    .groups = "drop"
  )

Joining with `by = join_by(month)`

# Calculate annual averages
annual_temps <- regional_clim %>% 
  mutate(n_days = days_in_month(date)) %>% 
  group_by(survey_area, year) %>% 
  summarise(
    sst = weighted.mean(sst, w = n_days),
    clim_sst = weighted.mean(clim_sst, w = n_days),
    sst_anom = weighted.mean(sst_anom, w = n_days),
    bot_temp = weighted.mean(bot_temp, w = n_days),
    clim_bt = weighted.mean(clim_bt, w = n_days),
    bt_anom = weighted.mean(bt_anom, w = n_days),
    .groups = "drop"
  )

Map of Areas

new_england <- ne_states("united states of america", returnclass = "sf")
canada <- ne_states("canada", returnclass = "sf")

(region_map <- regions %>% 
  mutate(area = fct_relevel(area, area_levels_long)) %>% 
  ggplot() +
  geom_sf(aes(fill = area), alpha = 0.6) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  theme_bw() +
  map_theme() +
  theme(legend.position = "inside", legend.position.inside = c(0.75, 0.18)) +
  scale_fill_gmri() +
  coord_sf(xlim = c(-76, -65.8), ylim = c(35.5, 44.4)) +
  labs(fill = ""))

Timeseries of Monthly Surface and Bottom Temperatures

The following plot shows what the monthly average temperatures have been based on these two datasets.

ggplot(regional_clim) +
  geom_line(aes(date, sst, color = "SST")) +
  geom_line(aes(date, bot_temp, color = "BT")) +
  scale_color_gmri(reverse = T) +
  facet_wrap(~survey_area, ncol = 1) +
  labs(y = "Temperature", color = "Color")

Warning: Removed 272 rows containing missing values or values outside the scale range
(`geom_line()`).

Warning: Removed 45 rows containing missing values or values outside the scale range
(`geom_line()`).

Timeseries of 1991-2020 Anomalies

Both surface and bottom span this 30-year period so we can view them as anomalies based on the same reference period.

.{panel-tabset}

Annual

ggplot(regional_clim) +
  geom_line(aes(date, sst_anom, color = "Surface Temperature Anomalies")) +
  geom_line(aes(date, bt_anom, color = "Bottom Temperature Anomalies")) +
  scale_color_gmri(reverse = T) +
  facet_wrap(~survey_area, ncol = 1) +
  labs(y = "Temperature Anomaly", color = "Color")

Warning: Removed 272 rows containing missing values or values outside the scale range
(`geom_line()`).

Warning: Removed 45 rows containing missing values or values outside the scale range
(`geom_line()`).

Seasonal

# Make these seasonal as if we were matching it to seasonal trawl data
# Calculate average for the season
seasonal_temps %>% 
  rename(`Bottom Temperature` = bt_anom, `Surface Temperature` = sst_anom) %>% 
  pivot_longer(names_to = "var", values_to = "vals", cols = ends_with("Temperature")) %>% 
  ggplot(aes(year, vals, color = survey_area)) +
  geom_point(alpha = 0.35, size = 0.5) +
  geom_ma(aes(linetype = "5-Year Moving Average"), n = 5, ma_fun = SMA) +
  facet_grid(fct_rev(var) ~ fct_rev(season)) +
  scale_color_gmri() +
  labs(y = "Temperature Anomaly\n(1991-2020 Baseline)",
       linetype = "",
       color = "Area")

Warning: Using the `size` aesthetic in this geom was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` in the `default_aes` field and elsewhere instead.

Warning: Removed 212 rows containing missing values or values outside the scale range
(`geom_point()`).

:::

Long-Term Trends and Net Change

These two datasets span a different range of years. For bottom temperature we will use from 1970 onward b/c we don’t have any other data before then in this project so we already are framing from 1970-onward anyways. Surface temperature will be done from 1982 onward because that is the first full year.

In the event of first order trends, net change over that time-span will be calculated as rate * years.

I’ll also be focusing on seasons to match the trawl survey as well.

Annual Trends

# Run the bottom tempoerature trends on annual summaries
start_yr <- 1970; end_yr <- 2019
annual_bottom_trends <- annual_temps %>% 
  filter(year %in% c(start_yr:end_yr)) %>% 
  split(.$survey_area) %>% 
  map_dfr(function(.x){
    annual_mod <- lm(bot_temp ~ year, data = .x)
    mod_summ <- summary(annual_mod)[["coefficients"]]
    
    # summary info
    mod_trend = ifelse(mod_summ[2,4] < 0.05, T, F)
    mod_rate = coef(annual_mod)[[2]]
    
    # tidy little output dataframe
    signif <- data.frame(
      "season" = c("Annual"),
      "trend" = c(mod_trend),
      "rate" = c(mod_rate)
    
      )
    return(signif)},
    .id = "survey_area") %>% 
  mutate(
    var = "Bottom Temperature\n(1970-2019)",
    temp_delta = ifelse(trend == T, rate * (end_yr - start_yr), NA))

# Run the bottom tempoerature trends on annual summaries
start_yr <- 1982; end_yr <- 2019
annual_surface_trends <- annual_temps %>% 
  filter(year %in% c(start_yr:end_yr)) %>% 
  split(.$survey_area) %>% 
  map_dfr(function(.x){
    # linear model
    annual_mod <- lm(sst ~ year, data = .x)
    annual_summ <- summary(annual_mod)[["coefficients"]]
  
    # summary info
    mod_trend = ifelse(annual_summ[2,4] < 0.05, T, F)
    mod_rate = coef(annual_mod)[[2]]
    
    # tidy little output dataframe
    signif <- data.frame(
      "season" = c("Annual"),
      "trend" = c(mod_trend),
      "rate" = c(mod_rate)
    
      )
    return(signif)},
    .id = "survey_area") %>% 
  mutate(
    var = "Surface Temperature\n(1982-2019)",
    temp_delta = ifelse(trend == T, rate * (end_yr - start_yr), NA))

# Combine those jokers
annual_trends_summary <- bind_rows(annual_surface_trends, annual_bottom_trends) %>% 
  mutate(survey_area = fct_relevel(survey_area, rev(area_levels_long)))

# Dumbell Plot
ggplot(annual_trends_summary) +
  geom_segment(aes(
    y = survey_area, yend = survey_area, 
    x = 0, xend = temp_delta, color = var), position = position_dodge()) +
  geom_point(aes(x = temp_delta, y = survey_area, color = var)) + 
  geom_vline(xintercept = 0) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  facet_wrap(~var, ncol = 2) +
  theme(panel.grid.major.y = element_blank()) +
  labs(
    x = expression(Delta~"Temperature"), 
    y = "", 
    title = "Net Temperature Change",
    fill = "")

Warning: Width not defined
ℹ Set with `position_dodge(width = ...)`

# Maybe this should be a table:

warming_table <- annual_trends_summary %>% 
  select(-c(season, trend)) %>% 
  rename(
    area = survey_area,
    `Net Change` = temp_delta, 
    `Rate of Warming` = rate) %>% 
  group_by(var) %>% 
  gt(rowname_col = "area") %>% 
  gt::tab_header(title = "Regional Surface and Bottom Temperature Change") %>% 
  fmt(columns = `Rate of Warming`, fns = function(x){str_c(signif(x,2), deg_c, " / year")}) %>% 
  fmt(columns = `Net Change`, fns = function(x){str_c(signif(x,2), deg_c)})
  

warming_table

Regional Surface and Bottom Temperature Change
	Rate of Warming	Net Change
Surface Temperature (1982-2019)
Gulf of Maine	0.044°C / year	1.6°C
Georges Bank	0.031°C / year	1.1°C
Southern New England	0.035°C / year	1.3°C
Mid-Atlantic Bight	0.02°C / year	0.74°C
Bottom Temperature (1970-2019)
Gulf of Maine	0.035°C / year	1.7°C
Georges Bank	0.024°C / year	1.2°C
Southern New England	0.02°C / year	0.97°C
Mid-Atlantic Bight	0.036°C / year	1.8°C

Seasonal (Spring & Fall) Trends

start_yr <- 1970; end_yr <- 2019
seas_bottom_trends <- seasonal_temps %>% 
  filter(year %in% c(start_yr:end_yr)) %>% 
  split(.$survey_area) %>% 
  map_dfr(function(.x){
    spring_mod <- lm(bot_temp ~ year, data = filter(.x, season == "Spring"))
    spring_summ <- summary(spring_mod)[["coefficients"]]
    fall_mod <- lm(bot_temp ~ year, data = filter(.x, season == "Fall"))
    fall_summ <- summary(fall_mod)[["coefficients"]]
    
    # summary info
    spring_trend = ifelse(spring_summ[2,4] < 0.05, T, F)
    spring_rate = coef(spring_mod)[[2]]
    fall_trend = ifelse(fall_summ[2,4] < 0.05, T,F)
    fall_rate = coef(fall_mod)[[2]]
    
    # tidy little output dataframe
    signif <- data.frame(
      "season" = c("Spring", "Fall"),
      "trend" = c(spring_trend, fall_trend),
      "rate" = c(spring_rate, fall_rate)
    
      )
    return(signif)},
    .id = "survey_area") %>% 
  mutate(
    var = "Bottom Temperature",
    temp_delta = ifelse(trend == T, rate * (end_yr - start_yr), NA))

# Run the bottom tempoerature trends
start_yr <- 1982; end_yr <- 2019
seas_surface_trends <- seasonal_temps %>% 
  filter(year %in% c(start_yr:end_yr)) %>% 
  split(.$survey_area) %>% 
  map_dfr(function(.x){
    spring_mod <- lm(sst ~ year, data = filter(.x, season == "Spring"))
    spring_summ <- summary(spring_mod)[["coefficients"]]
    fall_mod <- lm(sst ~ year, data = filter(.x, season == "Fall"))
    fall_summ <- summary(fall_mod)[["coefficients"]]
    
    # summary info
    spring_trend = ifelse(spring_summ[2,4] < 0.05, T, F)
    spring_rate = coef(spring_mod)[[2]]
    fall_trend = ifelse(fall_summ[2,4] < 0.05, T,F)
    fall_rate = coef(fall_mod)[[2]]
    
    # tidy little output dataframe
    signif <- data.frame(
      "season" = c("Spring", "Fall"),
      "trend" = c(spring_trend, fall_trend),
      "rate" = c(spring_rate, fall_rate)
    
      )
    return(signif)},
    .id = "survey_area") %>% 
  mutate(
    var = "Surface Temperature",
    temp_delta = ifelse(trend == T, rate * (end_yr - start_yr), NA))

# Combine those jokers
seas_trends_summary <- bind_rows(seas_surface_trends, seas_bottom_trends) %>% 
  mutate(survey_area = fct_relevel(survey_area, rev(area_levels_long)),
         season = fct_relevel(season, "Spring", "Fall"))

  
# Dumbell Plot
ggplot(seas_trends_summary) +
  geom_segment(aes(
    y = survey_area, yend = survey_area, 
    x = 0, xend = temp_delta, color = var), position = position_dodge()) +
  geom_point(aes(x = temp_delta, y = survey_area, color = var)) + 
  geom_vline(xintercept = 0) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  facet_grid(var~season) +
  theme(panel.grid.major.y = element_blank()) +
  labs(
    x = expression(Delta~"Temperature"), 
    y = "", 
    title = "Net Temperature Change",
    fill = "")

Warning: Width not defined
ℹ Set with `position_dodge(width = ...)`

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_segment()`).

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_segment()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_segment()`).

Warning: Removed 6 rows containing missing values or values outside the scale range
(`geom_point()`).

Baseline Shifting

This is more for conceptual figure. But to me, the broad-scale impacts to growth in a region are a manifestation of shifting baselines.

This could be shown with seasonal variability changes as a curve, or maybe as a heatmap.

Maybe just standardize all of them (monthly*regional averages) by the long term average?

# Scale on long-term average of each month
monthly_scaled_temps <- regional_temps %>% 
  group_by(survey_area, month) %>% 
  summarise(
    avg_sst = mean(sst, na.rm = T),
    avg_bt = mean(bot_temp, na.rm = T)) %>% 
  full_join(regional_temps) %>% 
  mutate(
    rel_sst = sst - avg_sst,
    rel_bt = bot_temp - avg_bt)

`summarise()` has grouped output by 'survey_area'. You can override using the
`.groups` argument.
Joining with `by = join_by(survey_area, month)`

Monthly Cycle Shifting Baselines

monthly_scaled_temps %>% 
  mutate(decade = floor_decade(year)) %>% 
  group_by(survey_area, decade, month = as.numeric(month)) %>% 
  summarise(
    rel_sst = mean(rel_sst, na.rm = T),
    rel_bt = mean(rel_bt, na.rm =T),
    obs_sst = mean(sst, na.rm = T),
    obs_bt = mean(bot_temp, na.rm = T)) %>% 
  pivot_longer(cols = starts_with("rel"), values_to = "rel_temp", names_to = "rel_var") %>% 
  pivot_longer(cols = starts_with("obs"), values_to = "obs_temp", names_to = "obs_var") %>% 
  mutate(
    rel_var = ifelse(str_detect(rel_var, "bt"), "Bottom Temperature", "Surface Temperature"),
    obs_var = ifelse(str_detect(obs_var, "bt"), "Bottom Temperature", "Surface Temperature"),
    month = str_pad(month, side = "left", pad = "0", width = 2),
    base_date = as.Date(str_c("2000-", month, "-01"))) %>% 
  ggplot() +
  geom_line(aes(base_date, obs_temp, color = decade)) +
  scale_color_gmri() +
  scale_x_date(date_breaks = "1 month", labels = date_format("%b")) +
  facet_grid(survey_area~obs_var)

`summarise()` has grouped output by 'survey_area', 'decade'. You can override
using the `.groups` argument.

Warning: Removed 72 rows containing missing values or values outside the scale range
(`geom_line()`).

Monthly Temperature Anomalies Heatmap

monthly_scaled_temps %>% 
  mutate(month = month.abb[as.numeric(month)],
         month = fct_rev(fct_relevel(month, month.abb))) %>% 
  select(-c(sst, bot_temp, avg_sst, avg_bt)) %>% 
  pivot_longer(names_to = "var", values_to = "rel_temp", cols = c(rel_bt, rel_sst)) %>% 
  mutate(var = if_else(var == "rel_bt", "Bottom Temperature", "Surface Temperature")) %>% 
  filter(var == "Bottom Temperature") %>% 
  drop_na(rel_temp) %>% 
  ggplot() +
  geom_tile(aes(x = year, y = month, fill = rel_temp)) +
  scale_fill_distiller(palette = "RdBu", limits = c(-4,4), oob = oob_squish, labels = label_number(suffix = deg_c)) +
  facet_grid(survey_area ~ var, scale = "free", space = "free") +
  guides(fill = guide_colorbar(
    barwidth = unit(5, "cm"), 
    frame.colour = "black",
    ticks.colour = "black")) +
  scale_x_continuous(expand = expansion(add = c(0,0))) +
  theme(legend.title.position = "top", 
        axis.text.y = element_text(size = 8),
  legend.title = element_text(hjust = 0.5)) +
  labs(fill = "Temperature Departure from All-Year Average",
       y = "")

Temperature Anomalies

Are regional anomalies (either annual or monthly) too correlated with actual temperatures to include together in a model?

# Get overall correlation
seasonal_temps %>% 
  select(bot_temp, bt_anom) %>% 
  drop_na() %>% 
  summarise(corr = cor(x = bot_temp, y = bt_anom))

# A tibble: 1 × 1
   corr
  <dbl>
1 0.291

seasonal_temps %>% 
  select(bot_temp, bt_anom) %>% 
  drop_na() %>% 
  ggplot(aes(bot_temp, bt_anom)) +
  geom_point() +
  labs(x = "Absolute Bottom Temperature", y = "Regional Bottom Temperature Anomaly", title = "Correlation: 0.29")

Modeling with Temperature & Anomalies

#### Load Data  ####
wigley_bmspectra_df <- read_csv(here::here("Data/model_ready/wigley_community_bmspectra_mod.csv")) 

Rows: 398 Columns: 17
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): season, survey_area, spectra_type, area_titles
dbl (13): est_year, xmin_fit, xmax_fit, xmin_actual, xmax_actual, n, b, conf...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Join temps Drop NA's
# Wrangle baby
wtb_model_df <-  seasonal_temps %>% 
  rename(est_year = year) %>% 
  mutate(survey_area = case_when(
    survey_area == "Gulf of Maine" ~ "GoM",
    survey_area == "Georges Bank" ~ "GB",
    survey_area == "Southern New England" ~ "SNE",
    survey_area == "Mid-Atlantic Bight" ~ "MAB")) %>% 
  left_join( 
    select(wigley_bmspectra_df, -bot_temp)) %>% #glimpse()
  drop_na(b, total_weight_lb, bot_temp, bt_anom) %>%
  group_by(survey_area, season) %>%
  mutate(
    roll_temp = zoo::rollapply(bot_temp, 5, mean, na.rm = T, align = "right",  fill = NA),
    roll_anom = zoo::rollapply(bt_anom, 5, mean, na.rm = T, align = "right",  fill = NA),
    .groups = "drop") %>% 
  mutate(
    yr_num = as.numeric(est_year),
    yr_fac = factor(est_year),
    survey_area = factor(survey_area, levels = area_levels),
    season = factor(season, levels = c("Spring", "Fall")),
    landings = total_weight_lb,
    yr_seas = str_c(season, est_year))

Joining with `by = join_by(survey_area, est_year, season)`

# What do they look like on the same x axis?
# Add median weight here:
wigley_medwt_df <- read_csv(here::here("Data/model_ready/wigley_community_medsize_mod.csv"))

Rows: 398 Columns: 16
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (3): survey_area, season, area_titles
dbl (13): est_year, n_species, numlen_adj_total, mean_len_cm, med_len_cm, mi...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

select(wigley_medwt_df, survey_area, season, est_year, med_len_cm, med_wt_kg) %>% 
  right_join(wtb_model_df) %>% 
  select(survey_area, season, est_year, med_len_cm, med_wt_kg, b, total_weight_lb, bot_temp, bt_anom) %>% 
  pivot_longer(
    cols = c(med_len_cm, med_wt_kg, b, total_weight_lb, bot_temp, bt_anom), 
    names_to = "var", values_to = "val") %>% 
  mutate(survey_area = fct_relevel(survey_area, area_levels)) %>% 
  ggplot(aes(est_year, val)) +
  geom_point(aes(color = survey_area), alpha = 0.35, size = 0.5) +
  geom_ma(aes(color = survey_area), n = 5, ma_fun = SMA, linetype = 1) +
  facet_grid(var~fct_rev(season), scales = "free") +
  scale_color_gmri() +
  theme(legend.position = "bottom", strip.text.y = element_text(angle = 0)) +
  labs(color = "5-Year Moving Average",
       y = "Value",
       x = "Year")

Joining with `by = join_by(survey_area, season, est_year)`

ggplot(wtb_model_df, aes(bt_anom, b)) +
  geom_point()

# Make some models using anomalies or temperatures

# actual mod
library(lmerTest)

Loading required package: lme4

Loading required package: Matrix

Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack

Attaching package: 'lme4'

The following object is masked from 'package:raster':

    getData

library(performance)
conflicted::conflicts_prefer(lmerTest::lmer)

[conflicted] Will prefer lmerTest::lmer over any other package.

# Lets compare which one does better
actual_temps_model <- lmer(
  b ~ survey_area * season * scale(roll_temp) + survey_area*scale(log(total_weight_lb)) + (1|est_year), 
  data = wtb_model_df)

anomalies_temps_model <- lmer(
  b ~ survey_area * season * scale(roll_anom) + survey_area*scale(log(total_weight_lb)) + (1|est_year), 
  data = wtb_model_df)

# # Make it make sense
# summary(two_temps_model)

check_model(actual_temps_model)

check_model(anomalies_temps_model)

performance(actual_temps_model)

# Indices of model performance

AIC      |     AICc |      BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma
--------------------------------------------------------------------------------
-812.907 | -809.940 | -727.170 |      0.670 |      0.665 | 0.016 | 0.062 | 0.064

performance(anomalies_temps_model)

# Indices of model performance

AIC      |     AICc |      BIC | R2 (cond.) | R2 (marg.) |   ICC |  RMSE | Sigma
--------------------------------------------------------------------------------
-788.602 | -785.634 | -702.864 |      0.670 |      0.665 | 0.016 | 0.062 | 0.064

# Quick plots
plot(ggeffects::ggpredict(actual_temps_model , terms = ~roll_temp*survey_area*season), add.data = T)

Warning: Argument `add.data` is deprecated and will be removed in the future.
  Please use `show_data` instead.

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

plot(ggeffects::ggpredict(anomalies_temps_model , terms = ~roll_anom*survey_area*season), add.data = T)

Warning: Argument `add.data` is deprecated and will be removed in the future.
  Please use `show_data` instead.

Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.